function [ population, EvaluationNum, GenerationNum ] = dsde(Problem, Algorithm)
%%
% PROBLEM SETTING
d = Problem.Dimension;
LB = Problem.LowerBound;
UB = Problem.UpperBound;
ObjFunc = Problem.ObjFunc;
% ALGORITHM ASSIGNMENT
PopulationSize = Algorithm.PopulationSize;
F = Algorithm.F;
CR = Algorithm.CR;
T = Algorithm.T;
N = Algorithm.MaxFEs;
%% INITIALIZATION
GenerationNum = 1;  % the current generation number
population = rand([PopulationSize,d]).*repmat((UB-LB),[PopulationSize,1])+repmat(LB,[PopulationSize,1]);
fitness = ObjFunc(population);
EvaluationNum = PopulationSize;  % the budget used
% max_fitness = max(fitness);
% min_fitness = min(fitness);
population = [fitness,population];  % [fitnesses, solution]
population_merge = zeros(2*PopulationSize, 1+d);
%% OPTIMIZATION
MutationVector = zeros(PopulationSize,d);
TrialVector = zeros(PopulationSize,d);
sc = ones(PopulationSize,1); % the number of iterations that the individual is not improved
while EvaluationNum<N
%% Mutation and Crossover
    % partiton the population into subpopulations
    GenerationNum = GenerationNum+1;
    population_merge(1:PopulationSize,:) = population;
    [ ClusterSize, ClusterNum, population ] = partition_population(population);
    id = 0;
    for i = 1:ClusterNum
        if i < ClusterNum
            ClusterSize_temp = ClusterSize;
        else
            ClusterSize_temp = PopulationSize - ClusterSize*(ClusterNum-1);
        end
        id_start = ClusterSize*(i-1)+1;
        % DE/lbest/1
        for j = 1:ceil(ClusterSize_temp/2)
            % mutation
%             r1 = ceil(rand().*ClusterSize_temp);
%             r2 = ceil(rand().*ClusterSize_temp);
%             while r1==r2
%                 r2 = ceil(rand().*ClusterSize_temp);
%             end
%             r1 = r1+id_start-1;
%             r2 = r2+id_start-1;
            r1 = ceil(rand().*PopulationSize);
            r2 = ceil(rand().*PopulationSize);
            while r1==r2
                r2 = ceil(rand().*PopulationSize);
            end
            id = id+1;
            MutationVector(id,:) = population(id_start,2:end) + F*(population(r1,2:end)-population(r2,2:end));
            % crossover
            u1 = (rand([1,d])<=CR);
            TrialVector(id,:) = population(id,2:end).*(1-u1) + MutationVector(id,:).*u1;
            u2 = ceil(rand().*d);
            TrialVector(id,u2) = MutationVector(id,u2);
            TrialVector(id,:) = min([TrialVector(id,:);UB]);
            TrialVector(id,:) = max([TrialVector(id,:);LB]);
        end
        % DE/current-to-rand/1
        for j = (ceil(ClusterSize_temp/2)+1):ClusterSize_temp
            % mutation and crossover
            r1 = ceil(rand().*PopulationSize);
            r2 = ceil(rand().*PopulationSize);
            while r1==r2
                r2 = ceil(rand().*PopulationSize);
            end
            r3 = ceil(rand().*PopulationSize);
            while r3==r1 || r3==r2
                r3 = ceil(rand().*PopulationSize);
            end
            id = id+1;
            TrialVector(id,:) = population(id,2:end) + rand([1,d]).*(population(r1,2:end)-population(id,2:end))...
                +F*(population(r2,2:end)-population(r3,2:end));
            TrialVector(id,:) = min([TrialVector(id,:);UB]);
            TrialVector(id,:) = max([TrialVector(id,:);LB]);
        end
    end
    % fitness
    population_merge((PopulationSize+1):end,2:end) = TrialVector;
    fitness = ObjFunc(TrialVector);
    population_merge((PopulationSize+1):end,1) = fitness;
    EvaluationNum = EvaluationNum+PopulationSize;
%% Adaptive Selection Mechanism
    [ clusters, seed ] = apc_cluster( population_merge );
    [ population, sc ] = selection(clusters, seed, population_merge, sc);
%% reinitialization
    restart_id = find(sc'>=T);
    if ~isempty(restart_id)
        for i=restart_id
            dis = sqrt(sum((repmat(population(i,2:end),[PopulationSize,1])-population(:,2:end)).^2,2));
            dis(i) = inf;
            [~,dis_id] = sort(dis);
            neighbor = dis_id(1:ClusterSize);
            worse_id = neighbor(population(neighbor,1)>=population(i,1));
            worse_id = [worse_id',i];
            NumWorse = length(worse_id);
            population(worse_id,2:end) = rand([NumWorse,d]).*repmat(UB-LB,[NumWorse,1])+repmat(LB,[NumWorse,1]);
            population(worse_id,1) = ObjFunc(population(worse_id,2:end));
            EvaluationNum = EvaluationNum+NumWorse;
            sc(worse_id) = 1;
        end
    end
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [ ClusterSize, ClusterNum, subpopulation ] = partition_population(population)
        % population: [subpopulation_id, fitness, solution]
        ClusterSize = floor(rand()*17+4);
        PopulationSize_f = size(population,1);
        ClusterNum = floor(PopulationSize_f/ClusterSize);
        distance = zeros(PopulationSize_f, PopulationSize_f);
        for ii = 1:PopulationSize_f
            distance(ii,:) = sqrt(sum((repmat(population(ii,2:end),[PopulationSize_f,1]) - population(:,2:end)).^2,2));
        end
        fitness_f = population(:,1);
        subpopulation = zeros(size(population));
        species_finish = 0;
        for ii = 1:ClusterNum
            [~, seed_id] = min(fitness_f);
            distance_temp = distance(seed_id,:);
            [~, sol_id] = sort(distance_temp);
            if ii == ClusterNum
                ClusterSize_f = PopulationSize_f - ClusterSize*(ClusterNum-1);
            else
                ClusterSize_f = ClusterSize;
            end
            species_id = (species_finish+1):(species_finish+ClusterSize_f);
            subpopulation(species_id,:) = population(sol_id(1:ClusterSize_f),:);
            [~,rank_id] = sort(subpopulation(species_id,1));
            subpopulation(species_id,:) = subpopulation(species_id(rank_id),:);
            fitness_f(sol_id(1:ClusterSize_f)) = inf;
            distance(sol_id(1:ClusterSize_f),:) = inf;
            distance(:,sol_id(1:ClusterSize_f)) = inf;
            species_finish = species_finish+ClusterSize_f;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [ cluster, seed ] = apc_cluster( population )
        lambda = 0.9;
        max_its = 100;
        con_its = 30;
        PopulationSize_f = size(population,1);
        distance = zeros(PopulationSize_f, PopulationSize_f);
        for ii = 1:PopulationSize_f
            distance(ii,:) = -sum((repmat(population(ii,2:end),[PopulationSize_f,1]) - population(:,2:end)).^2,2);
        end
        distance(1:(PopulationSize_f+1):end) = NaN;
        shared_value = median(distance,'all','omitnan');
%         shared_value = min(min(distance));
        distance(1:(PopulationSize_f+1):end) = shared_value;
        exemplar = [];
        a = zeros(PopulationSize_f,PopulationSize_f);
        r = zeros(PopulationSize_f,PopulationSize_f);
        its = 0;
        con_its_temp = 0;
        while its<max_its && con_its_temp<con_its
            its = its+1;
            con_its_temp = con_its_temp+1;
            r_old = r;
            a_old = a;
            % update resposibility
            as = a+distance;
            [as_max,as_max_id] = max(as,[],2);
            r_temp = repmat(as_max,[1,PopulationSize_f]);
            for ii = 1:PopulationSize_f
                r_temp(ii,as_max_id(ii)) = max(as(ii,[1:(as_max_id(ii)-1),(as_max_id(ii)+1):PopulationSize_f]));
            end
            r = distance-r_temp;
            % update availibility
            r_mod = r.*(r>0);
            sum_r_mod = sum(r_mod,1)-diag(r_mod)';
            a_temp = diag(r)'+sum_r_mod;
            a = repmat(a_temp,[PopulationSize_f,1])-r_mod;
            a = a.*(a<0);
            a(1:(PopulationSize_f+1):end) = sum_r_mod;
            % combine
            r = lambda*r_old + (1-lambda)*r;
            a = lambda*a_old + (1-lambda)*a;
            % construct clusters
            ar = a+r;
            exemplar_old = exemplar;
            [~,exemplar_i] = max(ar,[],2);
            exemplar = find(exemplar_i==((1:PopulationSize_f)'));
            if length(exemplar_old) ~= length(exemplar) || sum(abs(exemplar-exemplar_old))~=0
                con_its_temp = 0;
            end
        end
        % final results
        NumCluster = length(exemplar);
        seed = zeros(1,NumCluster);
        cluster = cell(1,NumCluster);
        for ii = 1:NumCluster
            solution_id = find(exemplar_i==exemplar(ii));
            [~,solution_id_rank] = sort(population(solution_id,1));
            cluster{ii} = solution_id(solution_id_rank);
            seed(ii) = cluster{ii}(1);
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [ new_population, new_sc ] = selection(clusters, seeds, population, sc)
%         {cluster} = [id]: from best to worst
        ClusterNum_f = length(clusters);
        seed_fval = population(seeds,1);
        P = (max(seed_fval)-seed_fval+1e-4)./(max(seed_fval)-min(seed_fval)+1e-4);  % probability
        new_population = zeros(PopulationSize, 1+d);
        new_sc = ones(PopulationSize, 1);
        
        new_population(1:ClusterNum_f,:) = population(seeds,:);
        new_sc_temp = find(seeds<=PopulationSize);
        new_sc(new_sc_temp) = sc(seeds(new_sc_temp))+1;
        id_f = ClusterNum_f;
        index = ones(ClusterNum_f,1).*2;
%         id_f = 0;
%         index = ones(ClusterNum_f,1);
        ii = 1;
        while id_f<PopulationSize
            if index(ii)<=length(clusters{ii})
                if rand()<P(ii)
                    id_f = id_f+1;
                    new_population_id = clusters{ii}(index(ii));
                    new_population(id_f,:) = population(new_population_id,:);
                    if new_population_id <= PopulationSize
                        new_sc(id_f) = sc(new_population_id)+1;
                    end
                    index(ii) = index(ii)+1;
                end
            end
            ii = ii+1;
            if ii > ClusterNum_f
                ii = 1;
            end
        end
    end
end

